/*
 * Grid.java
 *
 * Created on August 6, 2004, 10:42 PM
 */
package edu.umd.essic.maskbuilder;

import java.awt.*;
import java.awt.geom.*;
import java.lang.Math;
import java.util.Vector;
import java.io.*;

public class Grid {
    private float width, height;        //size in pixels
    private float step_x, step_y;       //step in pixels
    private int n_x, n_y;               //numbers of nodes
    private Point2D.Double[][] coord;    //array of koordinats
    private Point2D.Double[][] diff1;    // array of metrics
    private Point2D.Double[][] diff2;
    
   // variables for stretching
    private int stretch_type = 0;       // for stretch type traking
    private double[] xstr, ystr;         // resulted stretched coordinates 
    private float[] xreg, yreg;         // initial regular coordinates
    private float[] cen_p_x, cen_p_y, cen_x, cen_y;   // arrays of central point coordinates of stretched areas
    private float[] skew_x, skew_y, dim_x, dim_y;   // arrays of central point coordinates of stretched areas
    private int areas_num = 0;
 
    // misc. variables for stretching
    private float[] alpha_x, beta_x, gamma_x, alpha_y, beta_y, gamma_y;
    private boolean flag_equator = true;
    private double[] abc_x, abc_y;
    private float scalex, scaley;
    //private float skewx = 3;
    //private float skewy = 3;
    //private float dimx  = 30;
    //private float dimy  = 30;
    //private static float PI = 3.1415926535f;
    private double deltax, deltay;
    private float left_x, upper_y, nlon, nlat;
    private int Nmax;
    
    /** Creates a new instance of Grid */
    public Grid(int Width, int Height, int NX, int NY, float Left_x, float Upper_y, float Nlon, float Nlat) {
        width =  Width;
        height = Height;
        left_x = Left_x;
        upper_y= Upper_y;
        nlon = Nlon;// number of degrees in longitude
        nlat = Nlat;// number of degrees in lattitude
        n_x = NX;
        n_y = NY;
        coord = new Point2D.Double[n_x][n_y];
        diff1 = new Point2D.Double[n_x][n_y];
        diff2 = new Point2D.Double[n_x][n_y];
        for(int i=0;i<n_x;i++)
            for(int j=0;j<n_y;j++){
               diff1[i][j] = new Point2D.Double(1,1);
               diff2[i][j] = new Point2D.Double(0,0);
            }
        Build();
    }
    
    public void ChangeStep(int NX, int NY) {
        n_x = NX;
        n_y = NY;
        coord = new Point2D.Double[n_x][n_y];
        Build();
    }
    
    public void Build() {
        step_x = width / (n_x-1);
        step_y = height / (n_y-1);
        float right_x = left_x + nlon;
        float lower_y = upper_y- nlat;
        deltax = (right_x - left_x)/(n_x-1);  //longitudinal grid step in degrees
        deltay = (upper_y - lower_y)/(n_y-1); //lattitudinal grid step in degrees
        //coordinats X
        for(int ix=0;ix<n_x;ix++){
            for(int jy=0;jy<n_y;jy++){
               coord[ix][jy] = new Point2D.Double((left_x+ix*deltax),(upper_y-jy*deltay));
           }
        }    
     }
    public void Init2DStretching(int n){
      areas_num = n;
      cen_p_x = new float[areas_num];
      cen_p_y = new float[areas_num];  
      skew_x  = new float[areas_num];
      skew_y  = new float[areas_num]; 
      dim_x   = new float[areas_num];
      dim_y   = new float[areas_num];       
    }
    public void SetStretchedAreas(Point2D.Float[] c,Point2D.Float[] s, Point2D.Float[] d){
      if((c!=null) && (s!= null) && (d!=null)){ 
          for (int i=0; i<areas_num; i++){
                cen_p_x[i] = c[i].x; 
                cen_p_y[i] = c[i].y;
                skew_x[i]  = s[i].x;
                skew_y[i]  = s[i].y;
                dim_x[i]   = d[i].x;
                dim_y[i]   = d[i].y;                
          }
      }
    //  if(s!=null){ skewx = (float)s.getX(); skewy = (float)s.getY();}
    //if(d!=null){ dimx  = (float)d.getX(); dimy  = (float)d.getY();}
    }
    
    public void InitEssicStretching(int N){
      Nmax = N;
      alpha_x = new float[Nmax];  
      beta_x  = new float[Nmax];
      gamma_x = new float[Nmax];
      alpha_y = new float[Nmax];
      beta_y  = new float[Nmax];
      gamma_y = new float[Nmax];
    }
    public void SetEssicCoefficients(float[] AlphaX, float[] BetaX, float[]GammaX, float[] AlphaY, float[] BetaY, float[] GammaY, boolean flag){
       for(int i=0;i<Nmax;i++){ alpha_x[i] = AlphaX[i];}
       for(int i=0;i<Nmax;i++){ beta_x[i]  = BetaX[i]; }
       for(int i=0;i<Nmax;i++){ gamma_x[i] = GammaX[i];}
       for(int i=0;i<Nmax;i++){ alpha_y[i] = AlphaY[i];}
       for(int i=0;i<Nmax;i++){ beta_y[i]  = BetaY[i]; }
       for(int i=0;i<Nmax;i++){ gamma_y[i] = GammaY[i];}
       flag_equator = flag;
    }
    
    public void InitStretchVariables(){
        xreg = new float[n_x]; xstr = new double[n_x]; 
        yreg = new float[n_y]; ystr = new double[n_y]; 
        for(int i=0;i<n_x;i++){xreg[i]=(float)coord[i][0].getX();}
        for(int j=0;j<n_y;j++){yreg[j]=(float)coord[0][j].getY();} 
     }
    
    public void FinalizeStretchVariables(){
        cen_p_x = null; cen_p_y = null;
        skew_x  = null; skew_y  = null;
        dim_x   = null; dim_y   = null;        
        alpha_x = null; beta_x  = null; gamma_x = null;
        alpha_y = null; beta_y  = null; gamma_y = null;
        xreg = null;  xstr = null; 
        yreg = null;  ystr = null;
        System.gc();
    }
    
    public int GetEssicNmax(){
        return Nmax;
    }
    
    public int Get2dExpAreasNum(){
        return areas_num;
    }
    
    public void SaveEssicGridData(FileWriter out) throws IOException 
    {
        int flag = 0;
        if (flag_equator) flag = 1;
        out.write(String.valueOf(flag)); out.write("\t");
        out.write(String.valueOf(Nmax)); out.write("\t");
        out.write("\n");
        for (int i=0; i<Nmax; i++) {
            out.write(String.valueOf(alpha_x[i])); out.write("\t");
            out.write(String.valueOf(beta_x[i])); out.write("\t");
            out.write(String.valueOf(gamma_x[i])); out.write("\t");
            out.write(String.valueOf(alpha_y[i])); out.write("\t");
            out.write(String.valueOf(beta_y[i])); out.write("\t");
            out.write(String.valueOf(gamma_y[i])); out.write("\t");
            out.write("\n");
        }
    }
    
    public void SaveExp2dGridData(FileWriter out) throws IOException
    {
        out.write(String.valueOf(areas_num)); out.write("\n");
        for (int i=0; i<areas_num; i++){
            out.write(String.valueOf(cen_p_x[i])); out.write("\t");
            out.write(String.valueOf(skew_x[i])); out.write("\t");
            out.write(String.valueOf(dim_x[i])); out.write("\t");
            out.write(String.valueOf(cen_p_y[i])); out.write("\t");
            out.write(String.valueOf(skew_y[i])); out.write("\t");
            out.write(String.valueOf(dim_y[i])); out.write("\t");            
            out.write("\n");
        }
    }
    
    public void ReadEssicGridData(BufferedReader in) throws IOException
    { 
        String line = in.readLine();
        Pos p = new Pos();
        int flag = p.ReadInt(line);
        if (flag == 1) flag_equator = true; else flag_equator = false;
        Nmax = p.ReadInt(line);
        InitEssicStretching(Nmax);
        p.SetZero();
        for (int i=0; i<Nmax; i++) {
            line = in.readLine();
            alpha_x[i] = p.ReadFloat(line); 
            beta_x[i]  = p.ReadFloat(line);
            gamma_x[i] = p.ReadFloat(line);
            alpha_y[i] = p.ReadFloat(line);
            beta_y[i]  = p.ReadFloat(line);
            gamma_y[i] = p.ReadFloat(line);
            p.SetZero();
        }   
        p = null;
    }
    
    public void Read2dExpGridData(BufferedReader in) throws IOException
    {
        String line = in.readLine();
        areas_num = Integer.parseInt(line.substring(0,1));
        Init2DStretching(areas_num);
        Pos p = new Pos();
        for (int i=0; i<areas_num; i++) {
            line = in.readLine();
            cen_p_x[i] = p.ReadFloat(line); 
            skew_x[i]  = p.ReadFloat(line); 
            dim_x[i]   = p.ReadFloat(line);            
            cen_p_y[i] = p.ReadFloat(line); 
            skew_y[i]  = p.ReadFloat(line); 
            dim_y[i]   = p.ReadFloat(line);             
            p.SetZero();
        }
        p = null;
    }
    
    public void ReadGridData(BufferedReader in, int func) throws IOException
    {
        if (func == 1)
            ReadEssicGridData(in);
        else if (func == 2)
            Read2dExpGridData(in);
        else if (func == 3){
            ReadEssicGridData(in);
            Read2dExpGridData(in);
        }
    }
       
/////////////////////////////////////////////////////////
// Routine Stretch makes different types of stretching //
// choice funcnum as 0 - no stretching; 1 - 1D         //
// stretching with ESSIC function, 2,3 -               //
// two-dimensional exponential and combined functions. //
/////////////////////////////////////////////////////////
public void Stretch(int funcnum)
{
int i,j,k;
double sum;
double xx,yy;
double xmin,xmax,ymin,ymax;
  switch (funcnum)
  {
    case 1: //=== The ESSIC function 
    stretch_type = 1;
    xmin = coord[0][0].x;
    xmax = coord[n_x-1][n_y-1].x;
    abc_x = GG_Utilities.getESSIC_ABC(xmin,xmax,alpha_x,beta_x,gamma_x,Nmax,false);
    for (i=0; i<n_x; i++){
        xx = (xmin+i*(xmax-xmin)/(n_x-1));
        xstr[i] = GG_Utilities.Newton(xmin,xmax,xx,abc_x,alpha_x,beta_x,gamma_x,Nmax);
    }
    ymin = coord[0][0].y;
    ymax = coord[n_x-1][n_y-1].y;
    abc_y = GG_Utilities.getESSIC_ABC(ymin,ymax,alpha_y,beta_y,gamma_y,Nmax,flag_equator);
    for (j=0; j<n_y; j++){
        yy = (ymin+j*(ymax-ymin)/(n_y-1));
        ystr[j] = GG_Utilities.Newton(ymin,ymax,yy,abc_y,alpha_y,beta_y,gamma_y,Nmax);
    }
    
    for(i=0;i<n_x;i++){
        for(j=0;j<n_y;j++){
            coord[i][j].setLocation(xstr[i],coord[i][j].getY());
        }
    }
    for(j=0;j<n_y;j++){
        for(i=0;i<n_x;i++){
            coord[i][j].setLocation(coord[i][j].getX(),ystr[j]);
        }
    }
    
    /* // Control!!! 
    for(i=0;i<n_x;i++){
        for(j=0;j<n_y;j++){
            coord[i][j].setLocation(GG_Utilities.ESSICfunc(xstr[i],abc_x,alpha_x,beta_x,gamma_x,1),coord[i][j].getY());
        }
    }
    //GG_Utitities.getESSIC_ABC(ymin,ymax, alpha_y, beta_y, gamma_y,Nmax);
    for(j=0;j<n_y;j++){
        for(i=0;i<n_x;i++){
            coord[i][j].setLocation(coord[i][j].getX(),GG_Utilities.ESSICfunc(ystr[j],abc_y,alpha_y,beta_y,gamma_y,1));
        }
    }
    // End of control!!!*/
        
    
    break;
    ///---------------
      case 2: //=== 2D exponential stretching
        stretch_type = 2;
        xmin = coord[0][0].x;
        xmax = coord[n_x-1][n_y-1].x;
        ymin = coord[0][0].y;
        ymax = coord[n_x-1][n_y-1].y;
        coord[0][0].y = ymin;
        for (i=0; i<n_x; i++){
            xx = (float)(xmin+i*(xmax-xmin)/(n_x-1));
            for (j=0; j<n_y; j++){
                yy = (float)(ymin+j*(ymax-ymin)/(n_y-1));
                coord[i][j].x = GG_Utilities.Newton(xmin,xmax,xx,coord[i][j].y,cen_p_x,cen_p_y,skew_x,dim_x,dim_y,areas_num);
                coord[i][j].y = GG_Utilities.Newton(ymin,ymax,yy,coord[i][j].x,cen_p_y,cen_p_x,skew_y,dim_y,dim_x,areas_num);
            }
        }
        
        /* // Control!!!
         for (i=0; i<n_x; i++)
            for (j=0; j<n_y; j++)
                 coord[i][j].setLocation(GG_Utilities.EXPfunc(coord[i][j].x,cen_p_x,coord[i][j].y,cen_p_y,skewx,dimx,dimy,areas_num),
                                         GG_Utilities.EXPfunc(coord[i][j].y,cen_p_y,coord[i][j].x,cen_p_x,skewy,dimy,dimx,areas_num));*/
        break;
    ///--------------
      case 3: //=== Combibed 1D+2D stretching
        stretch_type = 3;
        xmin = coord[0][0].x;
        xmax = coord[n_x-1][n_y-1].x;
        ymin = coord[0][0].y;
        ymax = coord[n_x-1][n_y-1].y;
        abc_x = GG_Utilities.getESSIC_ABC(xmin,xmax,alpha_x,beta_x,gamma_x,Nmax,flag_equator);
        abc_y = GG_Utilities.getESSIC_ABC(ymin,ymax,alpha_y,beta_y,gamma_y,Nmax,flag_equator);
        cen_x = new float[areas_num];
        cen_y = new float[areas_num];
        for (int n=0; n<areas_num; n++){
            cen_x[n] = (float)GG_Utilities.ESSICfunc(cen_p_x[n], abc_x, alpha_x, beta_x, gamma_x, Nmax);
            cen_y[n] = (float)GG_Utilities.ESSICfunc(cen_p_y[n], abc_y, alpha_y, beta_y, gamma_y, Nmax);
        }
        for (i=0; i<n_x; i++){
            xx = (float) (xmin+i*(xmax-xmin)/(n_x-1));
            for (j=0; j<n_y; j++){
                yy = (float)(ymin+j*(ymax-ymin)/(n_y-1));
                coord[i][j].x = GG_Utilities.Newton(xmin,xmax,xx,coord[i][j].y,abc_x,alpha_x,beta_x,gamma_x,Nmax,skew_x,dim_x,dim_y,cen_x,cen_p_y,areas_num);
                coord[i][j].y = GG_Utilities.Newton(ymin,ymax,yy,coord[i][j].x,abc_y,alpha_y,beta_y,gamma_y,Nmax,skew_y,dim_y,dim_x,cen_y,cen_p_x,areas_num);
                //coord[i][j].x = GG_Utilities.Newton1(xmin,xmax,xx,coord[i][j].y,abc_x,alpha_x,beta_x,gamma_x,Nmax,skewx,dimx,dimy,cen_p_x,cen_p_y,areas_num);
                //coord[i][j].y = GG_Utilities.Newton1(ymin,ymax,yy,coord[i][j].x,abc_y,alpha_y,beta_y,gamma_y,Nmax,skewy,dimy,dimx,cen_p_y,cen_p_x,areas_num);
            }
        }
      break;
    }
  CalcMetrics();
}
private void CalcMetrics(){  
    switch (stretch_type)
    {
    case 1:
         for (int i=1; i<n_x; i++){
             xstr[i] = 0.5*(xstr[i-1]+xstr[i]);
             diff1[i][1].x = GG_Utilities.DiffESSIC(xstr[i], abc_x, alpha_x, beta_x, gamma_x, Nmax);
         }
         for (int j=1; j<n_y; j++){
             ystr[j] = 0.5*(ystr[j-1]+ystr[j]);
             diff1[1][j].y = GG_Utilities.DiffESSIC(ystr[j], abc_y, alpha_y, beta_y, gamma_y, Nmax);
         }
         for (int i=1; i<n_x; i++)
            for (int j=1; j<n_y; j++){
                diff1[i][j].x = diff1[i][1].x;
                diff1[i][j].y = diff1[1][j].y;
            }
        break;
    case 2:
        for (int i=1; i<n_x; i++)
            for (int j=1; j<n_y; j++){
                double x = 0.5*(coord[i-1][j].x+coord[i][j].x);
                double y = 0.5*(coord[i][j-1].y+coord[i][j].y);
                diff1[i][j].x = GG_Utilities.DiffExp(x, cen_p_x,y,cen_p_y,skew_x,dim_x,dim_y,areas_num);
                diff1[i][j].y = GG_Utilities.DiffExp(y, cen_p_y,x,cen_p_x,skew_y,dim_y,dim_x,areas_num);
            }
        break;
    case 3:
        // Calculate derivatives following the chain rule
        double fx,diffx,diff2x2;
        for (int i=1; i<n_x; i++)
            for (int j=1; j<n_y; j++){
                double x = 0.5*(coord[i-1][j].x+coord[i][j].x);
                double y = 0.5*(coord[i][j-1].y+coord[i][j].y);
                fx      = GG_Utilities.ESSICfunc(x, abc_x, alpha_x, beta_x, gamma_x, Nmax);
                diffx   = GG_Utilities.DiffESSIC(x, abc_x, alpha_x, beta_x, gamma_x, Nmax);
                diff1[i][j].x =  (GG_Utilities.DiffExp(fx, cen_p_x,y,cen_p_y,skew_x,dim_x,dim_y,areas_num)*diffx);
            }
        for (int i=1; i<n_x; i++)
            for (int j=1; j<n_y; j++){
                double x = 0.5*(coord[i-1][j].x+coord[i][j].x);
                double y = 0.5*(coord[i][j-1].y+coord[i][j].y);               
                fx      = GG_Utilities.ESSICfunc(y, abc_y, alpha_y, beta_y, gamma_y, Nmax);
                diffx   = GG_Utilities.DiffESSIC(y, abc_y, alpha_y, beta_y, gamma_y, Nmax);
                diff1[i][j].y =  (GG_Utilities.DiffExp(fx, cen_p_y,x,cen_p_x,skew_y,dim_y,dim_x,areas_num)*diffx);
            }
        break;
   }
}
/*
private void CalcMetrics(){  
    switch (stretch_type)
    {
    case 1:
         for (int i=1; i<n_x; i++){
             xstr[i] = 0.5*(xstr[i-1]+xstr[i]);
             diff1[i][1].x = GG_Utilities.DiffESSIC(xstr[i], abc_x, alpha_x, beta_x, gamma_x, Nmax);
             diff2[i][1].x = GG_Utilities.Diff2ESSIC(xstr[i],abc_x, alpha_x, beta_x, gamma_x, Nmax);
         }
         for (int j=1; j<n_y; j++){
             ystr[j] = 0.5*(ystr[j-1]+ystr[j]);
             diff1[1][j].y = GG_Utilities.DiffESSIC(ystr[j], abc_y, alpha_y, beta_y, gamma_y, Nmax);
             diff2[1][j].y = GG_Utilities.Diff2ESSIC(ystr[j],abc_y, alpha_y, beta_y, gamma_y, Nmax);
         }
         for (int i=1; i<n_x; i++)
            for (int j=1; j<n_y; j++){
                diff1[i][j].x = diff1[i][1].x;
		diff2[i][j].x = diff2[i][1].x;
                diff1[i][j].y = diff1[1][j].y;
                diff2[i][j].y = diff1[1][j].y;
            }
        break;
    case 2:
        for (int i=1; i<n_x; i++)
            for (int j=1; j<n_y; j++){
                double x = 0.5*(coord[i-1][j].x+coord[i][j].x);
                double y = 0.5*(coord[i][j-1].y+coord[i][j].y);
                diff1[i][j].x = GG_Utilities.DiffExp(x, cen_p_x,y,cen_p_y,skewx,dimx,dimy,areas_num);
		diff2[i][j].x = GG_Utilities.Diff2Exp(x,cen_p_x,y,cen_p_y,skewx,dimx,dimy,areas_num);
                diff1[i][j].y = GG_Utilities.DiffExp(y, cen_p_y,x,cen_p_x,skewy,dimy,dimx,areas_num);
		diff2[i][j].y = GG_Utilities.Diff2Exp(y,cen_p_y,x,cen_p_x,skewy,dimy,dimx,areas_num);
            }
        break;
    case 3:
        // Calculate derivatives following the chain rule
        double fx,diffx,diff2x2;
        for (int i=1; i<n_x; i++)
            for (int j=1; j<n_y; j++){
                double x = 0.5*(coord[i-1][j].x+coord[i][j].x);
                double y = 0.5*(coord[i][j-1].y+coord[i][j].y);
                fx      = GG_Utilities.ESSICfunc(x, abc_x, alpha_x, beta_x, gamma_x, Nmax);
                diffx   = GG_Utilities.DiffESSIC(x, abc_x, alpha_x, beta_x, gamma_x, Nmax);
                diff2x2 = GG_Utilities.Diff2ESSIC(x,abc_x, alpha_x, beta_x, gamma_x, Nmax);
                diff1[i][j].x =  (GG_Utilities.DiffExp(fx, cen_p_x,y,cen_p_y,skewx,dimx,dimy,areas_num)*diffx);
		diff2[i][j].x =  (GG_Utilities.Diff2Exp(fx,cen_p_x,y,cen_p_y,skewx,dimx,dimy,areas_num)*
                                Math.pow(diffx,2) + diff1[i][j].x*diff2x2);
            }
        for (int i=1; i<n_x; i++)
            for (int j=1; j<n_y; j++){
                double x = 0.5*(coord[i-1][j].x+coord[i][j].x);
                double y = 0.5*(coord[i][j-1].y+coord[i][j].y);               
                fx      = GG_Utilities.ESSICfunc(y, abc_y, alpha_y, beta_y, gamma_y, Nmax);
                diffx   = GG_Utilities.DiffESSIC(y, abc_y, alpha_y, beta_y, gamma_y, Nmax);
                diff2x2 = GG_Utilities.Diff2ESSIC(y,abc_y, alpha_y, beta_y, gamma_y, Nmax);
                diff1[i][j].y =  (GG_Utilities.DiffExp(fx, cen_p_y,x,cen_p_x,skewy,dimy,dimx,areas_num)*diffx);
		diff2[i][j].y =  (GG_Utilities.Diff2Exp(fx,cen_p_y,x,cen_p_x,skewy,dimy,dimx,areas_num)*
                                Math.pow(diffx,2) + diff1[i][j].y*diff2x2);
            }
        break;
   }
}*/

public void GetMetrics(double[][] mx1, double[][] mx2, double[][] my1, double[][] my2){
    for (int i=0; i<n_x; i++)
        for (int j=0; j<n_y; j++)
	{
            mx1[i][j] = diff1[i][j].x;
            mx2[i][j] = diff2[i][j].x;
            my1[i][j] = diff1[i][j].y;
            my2[i][j] = diff2[i][j].y;
        }
}
public void GetMetrics(double[][] mx, double[][] my){
    for (int i=0; i<n_x; i++){
        for (int j=0; j<n_y; j++){
            mx[i][j] = diff1[i][j].x;
            my[i][j] = diff1[i][j].y;
        }
    }
}
//result of this method - polygon from a grid containing the parameter point
/*    public Polygon GetPolygon(int x, int y){
        boolean flag=false;
        Polygon p = null;
        Point2D.Float[][] crd = Get_Im_Coords();
        for(int i=0;i<n_x-1;i++){
            for(int j=0;j<n_y-1;j++){
                p = new Polygon();
                p.addPoint(Math.round(crd[i][j].x), Math.round(crd[i][j].y));
                p.addPoint(Math.round(crd[i][j+1].x), Math.round(crd[i][j+1].y));
                p.addPoint(Math.round(crd[i+1][j+1].x), Math.round(crd[i+1][j+1].y));
                p.addPoint(Math.round(crd[i+1][j].x), Math.round(crd[i+1][j].y));
                flag = p.contains(x,y);
                if(flag) {break;}
            }
            if(flag) {break;}
       }
        return p;
    }
 */
   public Polygon GetPolygon(int i, int j, float z){
        Polygon p = null;
        p = new Polygon();
        p.addPoint(GetImX(i,j,z)+1,GetImY(i,j,z)+1);
        p.addPoint(GetImX(i,j+1,z)+1,GetImY(i,j+1,z));
        p.addPoint(GetImX(i+1,j+1,z),GetImY(i+1,j+1,z));
        p.addPoint(GetImX(i+1,j,z),GetImY(i+1,j,z)+1);
        return p;
    }
   public Polygon GetPolygon(int i, int j, float z, boolean grid_shown){
        Polygon p = null;
        p = new Polygon();
        if (grid_shown){
            p.addPoint(GetImX(i,j,z)+1,GetImY(i,j,z)+1);
            p.addPoint(GetImX(i,j+1,z)+1,GetImY(i,j+1,z));
            p.addPoint(GetImX(i+1,j+1,z),GetImY(i+1,j+1,z));
            p.addPoint(GetImX(i+1,j,z),GetImY(i+1,j,z)+1);
        } else {
            p.addPoint(GetImX(i,j,z),GetImY(i,j,z));
            p.addPoint(GetImX(i,j+1,z),GetImY(i,j+1,z));
            p.addPoint(GetImX(i+1,j+1,z),GetImY(i+1,j+1,z));
            p.addPoint(GetImX(i+1,j,z),GetImY(i+1,j,z));            
        }
        return p;
    }   
   public Polygon GetPolygon(int i, int j){
        Polygon p = null;
        p = new Polygon();
        p.addPoint(GetImX(i,j,1),GetImY(i,j,1));
        p.addPoint(GetImX(i,j+1,1),GetImY(i,j+1,1));
        p.addPoint(GetImX(i+1,j+1,1),GetImY(i+1,j+1,1));
        p.addPoint(GetImX(i+1,j,1),GetImY(i+1,j,1));
        return p;
    }
    public float getLonMin(){
        return (float)coord[0][0].x;
    }
    public float getLonMax(){
        return (float)coord[n_x-1][n_y-1].x;
    }
    public float getLatMin(){
        return (float)coord[0][0].y;
    }
    public float getLatMax(){
        return (float)coord[n_x-1][n_y-1].y;
    }
        
    public Point get_index(float longitude, float lattitude){
       Point index = new Point(0,0);
       switch (stretch_type){
           case 0:
               index.x = (int)((longitude-left_x)/deltax);
               index.y = (int)((upper_y-lattitude)/deltay);
               break;
           case 1:  
               index.x = (int)((GG_Utilities.ESSICfunc(longitude,abc_x,alpha_x, beta_x, gamma_x, Nmax)-left_x)/deltax);
               index.y = (int)((upper_y-GG_Utilities.ESSICfunc(lattitude,abc_y,alpha_y, beta_y, gamma_y, Nmax))/deltay);    
               break;
           case 2:  
               index.x = (int)((GG_Utilities.EXPfunc(longitude,cen_p_x,lattitude,cen_p_y,skew_x,dim_x,dim_y,areas_num)-left_x)/deltax);
               index.y = (int)((upper_y-GG_Utilities.EXPfunc(lattitude,cen_p_y,longitude,cen_p_x,skew_y,dim_y,dim_x,areas_num))/deltay);    
               break;    
           case 3:  
               index.x = (int)((GG_Utilities.EXPfunc(GG_Utilities.ESSICfunc(longitude,abc_x,alpha_x, beta_x, gamma_x, Nmax),cen_x,lattitude,cen_p_y,skew_x,dim_x,dim_y,areas_num)-left_x)/deltax);
               index.y = (int)((upper_y-GG_Utilities.EXPfunc(GG_Utilities.ESSICfunc(lattitude,abc_y,alpha_y, beta_y, gamma_y, Nmax),cen_y,longitude,cen_p_x,skew_y,dim_y,dim_x,areas_num))/deltay);    
               //index.x = (int)(GG_Utilities.ESSICfunc(GG_Utilities.EXPfunc(longitude,cen_p_x,lattitude,cen_p_y,skewx,dimx,dimy,areas_num),abc_x,alpha_x, beta_x, gamma_x, Nmax)-left_x);
               //index.y = (int)(upper_y-GG_Utilities.ESSICfunc(GG_Utilities.EXPfunc(lattitude,cen_p_y,longitude,cen_p_x,skewy,dimy,dimx,areas_num),abc_y,alpha_y, beta_y, gamma_y, Nmax));    
           
               break;
       }
       return index; 
    }
    
    public float xcoord(int i, int j){
        return (float)coord[i][j].x;
    }
    public float ycoord(int i, int j){
        return (float)coord[i][j].y;
    }    
   
    private int GetImX(int i, int j, float zoom) {
        double x = coord[i][j].x;
        return (int) Math.round(zoom*((x-left_x)/deltax)*step_x);
    }  
    
    private int GetImY(int i, int j, float zoom) {
        double y = coord[i][j].y;
        return (int) Math.round(zoom*((upper_y-y)/deltay)*step_y);
    }  
       
    /*public Point2D.Float GetNode(int i, int j){
        Point2D.Float crd = new Point2D.Float((float)(0.5*(GetImX(i,j)+GetImY(i+1,j))),(float)(0.5*(GetImY(i,j)+GetImY(i,j+1))));
        return crd;
   */
    
    // Get Graphical Coordinates 
    public Point2D.Float[][] Get_Im_Coords() {
        Point2D.Float[][] result = new Point2D.Float[n_x][n_y];
        for(int i=0;i<n_x;i++){
            for(int j=0;j<n_y;j++){
            result[i][j] = new Point2D.Float((float)(((coord[i][j].x-left_x)/deltax)*step_x), (float)(((upper_y-coord[i][j].y)/deltay)*step_y)) ;
            }
        }    
        return result;
    }    
    public Point[][] Get_Im_Coords(float zoom_level) {
        Point[][] result = new Point[n_x][n_y];
        for(int i=0;i<n_x;i++){
            for(int j=0;j<n_y;j++){
            result[i][j] = new Point((int)Math.round(zoom_level*((coord[i][j].x-left_x)/deltax)*step_x), 
                                     (int)Math.round(zoom_level*((upper_y-coord[i][j].y)/deltay)*step_y));
            }
        }    
        return result;
    }      
    // Get Actual Coordinates in degrees
    public Point2D.Double[][] Get_Re_Coords() {
        Point2D.Double[][] result = new Point2D.Double[n_x][n_y];
        for(int i=0;i<n_x;i++){
            for(int j=0;j<n_y;j++){
            result[i][j] = new Point2D.Double(coord[i][j].x,coord[i][j].y) ;
            }
        }    
        return result;
    }    
    public void SetCenteredGrid(){
       n_x--; n_y--;
        for (int i=0; i<n_x; i++){
            for (int j=0; j<n_y; j++){
                coord[i][j].setLocation(0.5*(coord[i][j].x+coord[i+1][j].x),0.5*(coord[i][j].y+coord[i][j+1].y));
            }
       }
    }

    public void SetNodeGrid(Point2D.Double[][] crd){
       n_x++; n_y++;
       for (int i=0; i<n_x; i++){
            for (int j=0; j<n_y; j++){
                coord[i][j].setLocation(crd[i][j].x,crd[i][j].y);
            }
       }
    }

    public float Get_StepX() {
        return (float)deltax;
    }  
    
    public float Get_StepY() {
        return (float)deltay;
    }  
    
    public int Get_NX() {
        return n_x;
    }    
       
    public int Get_NY() {
        return n_y;
    }    
    
    public float Get_Height() {
        return height;
    }    
    
    public float Get_Width() {
        return width;
    }    

    public int Get_StrType() {
        return stretch_type;
    }    
    
    public float Get_MinDX(){
        float[][] dx = new float[n_x][n_y];
        for (int i=1; i<n_x-1; i++)
	for (int j=1; j<n_y-1; j++){
		dx[i][j] = Math.abs((float)(coord[i][j].x-coord[i-1][j].x));
	}
        return minmass(dx);   
    }
    
    public float Get_MaxDX(){
        float[][] dx = new float[n_x][n_y];
        for (int i=1; i<n_x-1; i++)
	for (int j=1; j<n_y-1; j++){
		dx[i][j] = Math.abs((float)(coord[i][j].x-coord[i-1][j].x));
	}
        return maxmass(dx);   
    }
    
    public float Get_MinDY(){
        float[][] dy = new float[n_x][n_y];
        for (int i=1; i<n_x-1; i++)
	for (int j=1; j<n_y-1; j++){
		dy[i][j] = Math.abs((float)(coord[i][j].y-coord[i][j-1].y));
	}
        return minmass(dy);      
    }
    
    public float Get_MaxDY(){
        float[][] dy = new float[n_x][n_y];
        for (int i=1; i<n_x-1; i++)
	for (int j=1; j<n_y-1; j++){
		dy[i][j] = Math.abs((float)(coord[i][j].y-coord[i][j-1].y));
	}
        return maxmass(dy);      
    }
    
    private float minmass(float[][] mass){
        float mmm = 1.0e7f;
	for (int i=1; i<n_x-1; i++)
		for (int j=1; j<n_y-1; j++)
			if (mmm > mass[i][j]) mmm = mass[i][j];

        return mmm;
    }

    private float maxmass(float[][] mass){
        float mmm = 0.f;
	for (int i=1; i<n_x-1; i++)
		for (int j=1; j<n_y-1; j++)
			if (mmm < mass[i][j]) mmm = mass[i][j];

        return mmm;
    }

 

}